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pvi Abstract We study the long-term dynamics of a planetary system composed of 
, a star and a planet. Both bodies are considered as extended, non-spherical, rotat- 
ed ing objects. There are no assumptions made on the relative angles between the 
^^ orbital angular momentum and the spin vectors of the bodies. Thus, we analyze 
^^ full, spatial model of the planetary system. Both objects are assumed to be de- 
^~H formed due to their own rotations, as well as due to the mutual tidal interactions. 
^"^ The general relativity corrections are considered in terms of the post-Newtonian 
^_^ approximation. Besides the conservative contributions to the perturbing forces, 
n , there are also taken into account non-conservative effects, i.e., the dissipation of 
PtI the mechanical energy. This dissipation is a result of the tidal perturbation on the 
velocity field in the internal zones with non-zero turbulent viscosity (convective 
zones). Our main goal is to derive the equations of the orbital motion as well 
as the equations governing time-evolution of the spin vectors (angular velocities). 
O We derive the Lagrangian equations of the second kind for systems which do not 
.^ conserve the mechanical energy. Next, the equations of motion are averaged out 
^ over all fast angles with respect to time-scales characteristic for conservative per- 
I I turbations. The final equations of motion are then used to study the dynamics of 
the non-conservative model over time scales of the order of the age of the star. 
'"^ We analyze the final state of the system as a function of the initial conditions. 
J^ Equilibria states of the averaged system are finally discussed. 

'^1 Keywords Celestial mechanics • planetary system • energy dissipation 
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\_J 1 Introduction 

o 



Since the first discovery of an exoplanet revolving around the main sequence star 
(Mayor and Queloz|1995l, many planetary companions were detected in a few day 



orbits . The dynamics of such systems are strongly affected by relativistic as well 
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as non-point and non-Newtonian effects. These perturbations cause the periapses 
rotation and the precession of the orbital nodes as well as spins of the rotating 
bodies. Considering the longest time scale which is comparable to the age of the 
parent star, one has to take into account also a dissipation of the mechanical energy. 
It is believed, that the most important physical mechanism dissipating the energy 
is the tidal perturbation of the velocity field in parts of bodies possessing non-zero 
turbulent viscosity (e.g., 'Zahn"1977'). This takes place in the convective zones of 
stars and planets. The mechanisms of the energy dissipation, particularly active 



in stars and Jupiter-like planets, were studied by many authors (e.g., Goldreich 
and S oter"196 6| [Ogilvie and LinpOOl) |Wu||2005a|b| [Ogilvie and Lin|[2007| [Miller 



et al 2009: Gu and Ogilvie 



2009 



Arras and Socrates 20101. An open problem is 



to estimate values of physical parameters characterising the strength and time- 
scales of these processes, and in turn, the time-scale of the dissipative evolution of 
planetary systems. 

It is well known, that the energy loss leads to a variation (a decrease in general) 
both the semi-major axis and eccentricity, as well as to evolution of spins, their 
directions and magnitudes. The planetary dynamics of systems with energy loss 



were considered both in cases with one and two planets (e.g., Mardling and Lin 



2002i[Wltre and Savonije|2002|[Dobbs-Dixon et al|2004| |Mardlin g|2007[ [Barker aiid 
Ogilvie||2009l [Leconte et al||2010| [Rod riguez a nd Ferraz -Mello 2 0Tol [Michtchenko 



and Rodriguez|20il Rodriguez et al||2011[ iCorreia et al||2011[ [Laskar et al||2012l. 



The equations of motion of such systems are usually derived through introducing a 
dissipative force. For instance, the well known model of |Hut| ( [l981[ ), assumes that 
this force emerges due to a time delay in forming the tidal bulge and the orbital 
motion of a perturber. That implies a non-zero angle between the radius vector of 
the deformable body with respect to its companion, and the axis of symmetry of 
the tidally deformed object. In more elegant way, the tidal force was derived from 
the energy loss function E defined by Eggleton et al ( 1998 1 . 



In this paper, we found a more straightforward derivation of the dissipative 
model that relies on the Lagrangian equations of the second kind. As we will 
show, this approach makes it possible to obtain quite simply the dissipative forces 
acting in the A''-body system; however, we limit here the derivation to TV = 2. 
Moreover, to derive the equations governing the evolution of angular velocities, 
both in conservative, as well as in non-conservative models, we should not apply 
the Euler equations, which hold only for a specific form of the potential energy V 
that has to be then a function of Euler angles (p, 9, tp, and should not depend on 
their time derivatives 0, 9, ijj- This is only true in the case of the rigid body. For 
deformable objects, V is a function of the angular velocity O = 0{(f>, 9, ip, (j), 9, ip) 
and thus it does not fulfill these assumptions. Hence, the Euler equations stating 
that the time derivative of the rotational angular momentum equals to the torque 
acting on the rigid body do not hold in general. However, we will show here that 
it is still possible to obtain the equations of the evolution of the angular velocities 
in vectorial form, which is reminiscent of the classic Eulerian equations. 

The plan of this paper is the following one. In Section 2 we derive the equations 
of motion for a general form of the Lagrangian C = To{Oo) +Ti{ni) — Vo(r,i?o) — 
yi(r, Qi)+£,i{v,r), and a dissipative function E = Eo{r, r, no) + Ei{r, r,/2i). Here, 
symbols r, r denote the planetary position and the orbital velocity vectors relative 
to the star, and fioi'f?! stand for the angular velocity vectors of the star and the 
planet, respectively. 
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In the next Section 3, we find expressions for these functions in a particular 
model considered in this paper. We use the polytropic model of |Chandrasekhar| 
( 1933a|b ) to calculate the internal structure and a deformation of figures of both 
objects. Relativistic correction to the Lagrangian is taken from Brumberg (20071. 



To determine the energy loss function, we use a simple model by Eggleton et al| 
(Il998l). 



In Section 4, the derived equations of motion are averaged out over all angles 
that vary in time-scales related to the conservative evolution. These angles, or- 
dered from the fastest to the slowest one are the following: the mean anomaly, 
the horizontal angle of the precessing planetary spin in the orbital frame, and the 
argument of pericenter. As the result we obtain the equations of motion describing 
the dissipative dynamics only, which are then studied in Section 5. We analyze the 
final state of the planetary system in terms of the initial conditions. As a particular 
solution, we discuss the equilibrium permitted in the system. 



2 The equations of motion 

We shall consider a planetary system in terms of the mechanical system with the 
Lagrangian £ = C{qi, . . . ,qn,qi, ■ ■ ■ ,qn) and energy loss function E = E{qi, . . . ,qn,qi, 
where q^ are the generalized coordinates, and g,i are the generalized velocities. In- 
dex i spans 1 to n, where n is the number of the degrees of freedom. The dynamical 
evolution of the system is governed by solutions to the Lagrangian equations of 



the second kind (e.g., Greiner (2003), p. 328) 



d 
dt 



dC 
dqi 



dC 
dqi 



IdE 
2d4i' 



. ,n. 



(1) 



The above equations are correct when E fulfills the following condition (e.g.. Gold- 
stein et al||2002[ p. 63) 

j:^^IS.=^e. (2) 



dqi 



Particularly, it takes place when _B is a quadratic form of the generalized velocities. 
The explicit form of E is given further in the text. Although it is not a quadratic 
form in q, it fulfills the above condition as well (it is shown in Appendix [A| . 

In the problem considered here, the generalized coordinates are the following: 

{<li}l=l = {x,y,Z,Sx,Sy,Sz,Px,Py,Pz}- 

The first three coordinates are components of the position vector of the planet 
relative to the star, r = [x,y,z\ . The last six coordinates are vectors of three in- 
dependent components in two unit quaternions, s = [sx,Sy, Sz] , p = [px,Py,Pz\ '■ 

3 = So+isi+jS2+tS3, Sfe G R, 
p=P0+ipi+iP2+ip3, Pfc G R- 

The components of the quaternions as well as their time derivatives are related as 
follows: 

So + Sl + S2 + S3 = 1, So So + Si Si + S2 S2 + S3 S3 = (3) 



,<?"), 
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and similarly, for p. Therefore, only three components of each quantity are indepen- 
dent. We choose these independent components as follows: s^ = si,Sy = S2,Sz = S3 
and px = pi,Py = P2,Pz = P3- Thus so, so,po,Po are functions of s,s, p,p, according 
to Eq. ([3|. 

The Lagrangian of the system, as well as E, are assumed to be functions of plan- 
etary position, the velocity and spin vectors of both bodies, i.e., C = C{r,r,no,^i), 
E = E(v,v,Oq,Qi). It is well known, that the angular velocity may be expressed 
with the help of quaternions in the following form, e.g.. Heard (2006), p. 49: 







/ So Si 


S2 


S3 \ 


no) 


= 2 


—Si so 

-S2 S3 


-S3 

So 


S2 

-Si 






\-S3 -S2 


Si 


SO / 



fso\ 

Si 

S2 

V^3/ 



(4) 



and similarly for Oi with quaternion p. Thus Oq 
There exists also an inverse relation: 



l2o(s,s) and Oi = /2i(p,p). 



fso\ 

Sl 

S2 

VS3/ 



f So -Sl -S2 -S3^ 

51 So S3 — S2 

52 -S3 So Sl 
\S3 S2 — Sl So / 





Oo 



(5) 



The angular velocities are then 17o = ■f2o(s, s) and Oi = i?i(p,p). 
The full set of Lagrange equations are then the following: 



d 
dt 
d 
dt 
d_ 
dt 



d£ 
dr 
dC 
ds 
dC 
dp 



dC _\ dE 
dv 2 dr ' 
d_C_ldE_ 
~ds ~ 2 ~dl' 
dC _l dE 
9p ~ 2 di>' 



(6) 



It may be shown (see AppendixlBl, that the second equation in the above set may 
be transformed into the matrix equation 



iX 



x = 



d 

dt 



dC 
dQo 



1 dE 

2d no' 



Y = 



dC 

dno' 



(7) 



2 I 2 

so + si 

-S0S3 -I- S1S2 

S0S2 + S1S3 



S0S3 + S1S2 

2,2 
Sq + S2 

-SoSl + S2S3 



-SqS2 + S1S3 

SlSQ + S2S3 

„2 I „2 
So + S3 



(8) 



(SoSo + SlSl 
— S0S3 -I- S2S1 
S0S2 + S3S1 



S0S3 + S1S2 — S0S2 + S1S3 

SOSO + S2S2 SoSl -|- S2S3 
— SoSl + S3S2 SoSo + S3S3 



(9) 



This equation may be solved with respect to X, leading to the following solution: 



X: 



(a^i A2) 



—^0,2 ^0,y \ 

f^o.z —f^o.x Y. 

— f^O,y ^0,x y 



(10) 
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where the angular velocity vectors have components Oq = [fio^x, ^o,v, ^o,z] for 
the star. Similar expression may be obtained for the angular velocity of the planet, 
fii = [Oi^x, ^i,y, ^i,z] ■ The product in the right-hand side of the equation is 
simply the vector product Qq x Y. The final equations of the evolution of /2; have 
the following form: 

The above equations are valid for systems, which Lagrangian depends on the 
space orientation of extended objects only through the angular velocities. As we 
will show, for the case considered here, these equations have similar explicit form 
as the Euler equations. Nevertheless, they were derived under assumption of a 
particular form of the potential function V = V{r, Oq, Oi), suitable to our model. 



3 The Lagrangian and the dissipative function 

The Lagrangian of the system is given by the following expression: 

£ = Tp-p - \/p-p - Vfot - Hid + ^rot + -Creb (12) 

where 

J^p-p-^P^' Vp U'^J 

are the kinetic and potential energies of two point masses interacting with accord 
to the Newtonian gravity. Masses of the star and the planet are denoted with mp 
and mi, respectively, /3 = (1/mo + l/m\)~^ is the reduced mass, r = ||r|| and k is 
the Gauss constant. Terms V^ot s-nd V^^i,^ are for the perturbing potential energy of 
two extended, non-spherical objects. We assume, that each object is deformed due 
to its own rotation, as well as to the mutual tidal interaction with the other body. 
These two effects are considered separately. Such a simplification is correct to the 
first order. Moreover, we assume that the planet deforms the shape of the star due 
to the point mass interaction, and the star is a point-mass perturber deforming 
the figure of the planet. The rotational kinetic energy of a deformable object is 
given by Trot) while the relativistic term of the Lagrangian is denoted by -Cj-gl- 

3.1 Potential of the axially symmetric object 

Both perturbing terms V^ot ^i^d V^j^j have an axial symmetry. The axis of symmetry 
of V^ot coincides with a direction of the angular velocity vector of a particular 
object, while the axis of symmetry of Vj;^ coincides with a direction of a vector r 
joining the mass centers of the bodies. It is well known, that the potential energy of 
a system that consists of a point mass m and extended axially symmetric object 
of mass M and characteristic radius R may be developed in Taylor series with 
respect to (R/r): 

^axial = ^^E'^2i - P2dv-z)^^^J2R^P2{r-i), (14) 
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where J21 are the Stokes coefficients, and P21 {x) are the Legendre polynomials [see, 
Schaub and Junkins] p003t , p. 480 or [Murray and Dermott] ( |2000[ ), p. 133]. 



e.g. 

Here, these series are truncated at terms proportional to J2. In the case of the 
rotational deformation, z = O/Q, where fi = \\0\\. For the tidal perturbations 
induced by a body at r, z = f = r/r. The zonal Stokes coefficients J2 are expressed 



by the following integrals (e.g., Schaub and Junkins||2003[ p. 477): 



rR'W 



J2 = -— ^ / / p(7^,■!9)7^^ sin1?P2(cos19)d19d7^ 

™^ ./o Jo 



.7 ./ 

TT i-R 

p{Tl,'d)n'' ain-d P2 {cos {f)Mdn, (15) 






.7 ./ 



where the density p {TZ, ■&) depends on the magnitude of perturbations due to the 
rotation and tides. That function may be determined in a coordinate system which 
has its 2-axis along z. Angle d is measured from this axis, between (the "north 
pole") to TT (the "south pole"), and the radial variable is 7?. G [0, 7?'(i9)]. The shape 
of a body is described in terms of -R'(i?), which may be approximated by R when 
the perturbation is small. 

3.1.1 Calculating the Love numbers 

To determine and describe tidal perturbations, we shall need to ca lculate the Love 
numbe rJM To accomplish this task, we apply remarkable results of Chandrasekhar 



(1933ap, who considered uniformly rotating polytropeal and obtained equations 



for the density function p(7?., i9). He postulated the following form of p: 

p(z,i3) = pcW'iz,^), z = ATZ, (16) 

where pc is the central density and n is the polytropic index and ^ is a function. 



which is usually introduced to define dimensionless variable z (see, e.g., Kippen- 
hahn and Weigert|1994 p. 176 for the explicit formulae). He found that for slowly 



rotating body: 

W {z, ^) = w{z) + -—2— [Ho (z) - an U2 (2) P2 (cos d)] . (17) 

Function w{z) may be obtained by solving the Lane-Emden equation (Eq. [12] in 
the cited paper). Functions Uo {z) and U2 {z) may be found by solving equations 
[37i] and [37 2] in the cited paper, respectively (let us note that the notation 
used here is different from the original one). Coefficients an depend on the poly- 
tropic index n and may be derived as solutions to the problem considered by 



^ In stellar astrophysics, the so called apsidal motion constant (instead of the Love number, 
which is twice as much) is used as a physical parameter giving the rate of the rotation of 
apsidal line of the binary orbit in space. On the other hand, in planetary astrophysics, the 
Love number is usually used (not only when a planet-moon system is considered, but also for 
a star-planet system). In this work we use this planetary convention. 

^ A polytrope is understood here as a gaseous object being in hydrostatic equilibrium under 
its own gravity. The pressure-density relation is given by the formulae P = ii'pt^+l)/" (where 
pressure is denoted by P, density by p, ii" is a constant factor and n is a polytropic index). 
A rotating polytrope or/and being attracted by some external force is then called a perturbed 
polytrope. 
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Chandrasekhar. Using the resulting p{TZ,i9), it is relatively easy to show that the 
Stokes coefficient J2 defined by Eq. ( 15 ) has the following form 



J2 = 

'^L.r 



,(rot) 



1 R^ n"^ 

3 L>l-^fc2^' 



«L,r 



(n) 



6a„ 

54 



n [w (z)] U2 {z) z dz, 



(18) 
(19) 



where Zn is dimensionless size of undisturbed body, i.e., w {zn) = and kj^ j. may 
be attributed to the fluid Love number of the object deformed due to the rotation 



(Munk and MacDonald ( 1975 I, p. 26), which is a function of index n. In the range 
n G [0,4], this function is very well approximated by the formulae: 



logio k 



10 '=L,r ~ .f ("-) = logio 7, +ain + a2n + 03 n 



Ql 



-0.4872, Q2 = +0.0424, 03 



-0.0238. 



(20) 



In the next paper, Chandrasekhar (1933b) considered the deformation of the poly- 



tropic body having mass M by the tidal force emerged due to the outer point-mass 
perturber of mass m. He found the density function, which may be then used to 
calculate J2 coefficient of tidally distorted object. It may be shown that for small 
perturbations, the linear approximation is valid, and then: 



J , (tidal) 

J2 - Jo 



M 



«L.t> 



(21) 



where fcn is the tidal Love number (Munk and MacDonald (19751, p. 27) and is 
given by the formulae: 



'i^L.t = fcL,t("-) 



5 Zn 

6 an 



[3W2 {zn) + ZnU2 ' {zn) 



L.r 



.(«) 



(22) 



Thus the fluid Love numbers of rotationally and tidally deformed object are nu- 
merically equal. Let's denote them by ki^. It is not surprising, because k\^ is a 
physical measure of deformability of an object under the attraction of some force 
(here, centrifugal and tidal forces are considered). The Love number relates to the 
quantity Q introduced in (Eggleton et al (1998), Eq. 15c), i.e., fc^ = Q/(l — Q). 
The numerical results obtained above are in agreement with the results stated in 
the cited work. However, the computation of fc^ (or the apsidal motion constant) 
were performed by many authors before (e.g., Brooker and Olle 1955), to make 
this work self-contained, we present a brief overview of one of the way which leads 
to numerical values of ki^. We choose to make use of Chandrasekhar's results, nev- 



ertheless this is not the only approach possible (see, e.g., Eggleton et al 19981. 



Equation 15c from this last paper gives a formulae for Q for a more general model 
of mass distribution. For a polytropic model they find a polynomial expression, 
Eq. 19. Our numerical results arisen from Eq. (20) are in agreement with the 
results in the literature (e.g., Brooker and 011e||1955 



Eggleton et al 19981 



However, the model of a polytrope being deformed by centrifugal as well as 
tidal forces is expected to work well for stars and gaseous planets, it is not expected 
so, when considering rocky planets (see, e.g., Bursa|[l984 ) . 
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To conclude, we write down the following forms of Vfot s-nd V^j,-! 

mifcLoi?gr2^ . mo fcL_i i?f r2? , 

14ot = 5-^ P2 r-r2o + ^-^ P2 r-/2i), (23 

fc^ m? 7?^ fcy n fc^ mn -R? fci i 

Hid = V^^ V^^' (24) 

where kj^ q and /cl 1 are the Love numbers of the star and the planet, respectively. 



Their numeric values are given by Eq. (20). It is believed, that a mass distribution 



in Sun- like stars is well approximated by a polytropic model of index n G [3,4]. In 
a case of Jupiter-like planets, the appropriate value of n £ [1,2]. 



3.2 The kinetic energy of rotating, extended bodies 

The kinetic energy of rotations of extended objects is given by a sum: 

T,ot = lo^hOo + ^nUiOi, (25) 

where Ip and Ii are the moments of inertia of the star and the planet, respectively. 
The moment of inertia is defined in a coordinate system with the origin that 



coincides with the mass center of the body as the following integral (e.g., Schaub 



and Junkins| ( |2003[ ), p. 129): 

ll= f piiv){-[r][r])d% (26) 

Jbody 1 

where f points towards the mass element in the body, and [f] is the so called tilde 
matrix of a vector f . For a vector a = [ax, ay,az] , it has the form of: 




(27) 

-i,y ax / 

The matrix multiplication of [a] and a vector b is equivalent to the the vector 



product a X b. A multiplication in Eq. ( 26 ) gives 



-„'2 I ;;2 



- [f] [r] = -xy i^ + -z^ ~y~z . (28) 

y —xz ^y z X +y J 

The density function of the l-th object {I = 0, 1) may be written as a sum: 

Pl{r)=p[°\r) + Apf'''\r), (29) 

where pj (f) is for the density function of undisturbed body, f = ||r||, ZipJ (?) 
is a correction stemming from the rotational deformation. Thus the moment of 
inertia may be expressed as the following sum: 

Ii = l|°^ -F Z\lf °*\ (30) 



K. = K(n) = 


2 2 

II = -^miRi Kn, 


34|™'(.n)l./o " ^"^" "" 
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Each term in that equation represents the moment of inertia for some density 
function. The first term is the moment of inertia of undistorted body possessing 
spherically symmetric distribution of its mass. Using a simple polytropic model of 
that object, we obtain: 

(31) 
(32) 

where E is a unit 3x3 matrix, i.e., E = diag(l, 1, 1). For the homogeneous mass 
distribution (n = 0), coefficient k = 1. In the range of n e [0,4], the coefficient 
is well approximated by the formulae similar to that ones derived for the Love 
number: 

logiQK ^ /(n) = Qi n + Q2n^ +a3'i^, (33) 

Qi = -0.2061, Q2 = +0.0297, as = -0.0140. (34) 

The rotational contribution to the moment of inertia in the principal axes frame 
(with the z axis determined by Q), has the following form: 

^f(rot) ^ ^(rot,l) jg _ ^(rot,2) ^.^^ ^^^ ^^ _^^^ ^3^^ 






(rot,l) ^ 4 nf Rf Ti (rot,2) ^ ^l Rj K,l 

' 3fc2 ' ' 9fc2 ■ 



(36) 



The coefficient r is given in the polytropic model by the integral: 

r = r{n) = \ f " nw'''\z)Uoiz) z* dz. (37) 

^n Jo 

Similarly to the coefficients k-^ and k, r may be also approximated by the polyno- 
mial 

log^QT Ri /(n) = Qi n + 02 71^ + as"-^, (38) 

Qi = -0.3818, Q2 = +0.0220, as = -0.0125. (39) 

The resulting form of Tj-ot is then: 

Trot = ^ (^0 + C) ^0 + i (/i + in) n\- (40) 

3.3 The relativistic perturbing Lagrangian 

The model considered here includes also the relativistic perturbations. We do not 
include terms containing Qq,Q\, because their contributions to the equations of 
motion are of the second order. The relativistic contribution to the Lagrangian of 
the two point masses is given by the formulae (e.g., Brumberg||2007 1 : 



1 f /•2\2 r^ (r-r)^ 1 1 

^rel=^|7l(r'] + ^2- + 73^^ - 74^ |, (41) 



10 
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where the mass parameters have the foUowing form 

_ /3 mg + ml k 

8 (mo + mi) 

73 = -ttP , 74 



72 = -^ f 3 mo mi +;8^j 



1« 2 



(42) 
(43) 



3.4 The dissipative function 

A function describing dissipation of the mechanical energy of the system is the 
following sum: 

E = Eo + Ei, (44) 

where the first term comes from the energy dissipation in the star due to tidal 
interaction with a planet, while the second term expresses the energy dissipation 
in the planet due to the interaction with the stax. The particular term of E are 



given by a simple model proposed by (Eggleton et alp 998 Eq. 43) 



E, = - 



9aimfRi^kl, 
2r8 



(r- 



+ r^ - 2 r • (fi; X r) + (/?, X r)' 



(45) 



for I = 0,1. Parameters ao, ai are dissipation constants for th e star and the planet, 
respectiveljrj They may be expressed in the following form (Kiseleva et al_^1998[): 



o"/ 



miRf 



1 + fcL,, 



k"^ mi 



(46) 



where Ao and Ai are non-dimensional coefficients of the energy loss rates in the 
star and the planet, respectively. A calculation of these values is complex, and 
we postpone the derivation to other work, actually, that is basically beyond the 
scope of the present paper. Following the literature (e.g., Ogilvie and Lin|2004 Wu 



2005a|b Ogilvie and Linj 2007| Hansen 2010), we may assume these coefficients 



in the range of Ao € [10~', 10"*] and Ai G [10~ , 10~ ]. Particular values of the 
dissipative coefficients are not well known. In the paper by Hansen ( |2010[ ), he 
calibrated the equilibrium tide model comparing the results of calculations with 
the observed relation of the orbital period and eccentricity of close binaries for 
which their age is known. He obtained a* « 5.3 x 10~^, which may be "translated" 
to Ao, through relation Ao = (t*(1 + fcLo)^^' which then gives Ao ~ 5 x 10^^. 
The dissipation constants for Jupiter-like planets were derived by analysis the 
period-eccentricity relations of known systems with these planets. The resulting 
o-p ^ 6.8 X lO^'^ gives Ai ^ 5 x lO^'^. 



^ We assume that ct; are constant, which is physically wrong, while these quantities depend 
on the tidal frequency (which is a function of orbital elements and time). Nevertheless, because 
of poor knowledge of their values, we m ake t his simplification of the model. For an overvie w on 
the more realisti c tidal models see, e.g., ( |Fer raz-Mello et al 2008 Efroimsky and Williams 2 009| 
|Efroimsky|20 Ilk . The dependence of cri on the tidal frequency is also discussed in subsection 4.3. 
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4 The explicit form of the equations of motion 

Furthermore, we consider a particular form of the Lagrangian and the dissipative 
function: 

L = To{no)+Ti{ni) - l/o(r,/2o) - Vi(r,fii) + ri(r, r), (47) 

E = Eo{r,r,nQ) + Ei{r,r,ni). (48) 

The derivative of £. over r does not depend on fii and similarly, the derivative of £ 
over particular angular velocity does not depend on the other n, nor on the linear 
velocity. Hence, the vector equations for r, Hq and i?i are separable. Equations in 



the top row of ([6| describe the orbital motion of the planet. The equations (111 
with I = 0,1 are for the evolution of the angular velocities of the star and the 
planet, respectively. 



4.1 The equations of translational motion 

The top-row vector equation of ([6| has to be solved with respect to r to obtain the 
explicit form of the equations of translational motion. The dependence of C on r 
appears in the Tp-p and £j.gj terms. Only the first one is a uniform quadratic func- 
tion of r. Nevertheless, we may solve the problem. The solution has the following 
form in the first order approximation: 

fi 6ii f Ao ^ Ai\ 
r ~ =-r =- \ r 



1 1 

'2P 



iE^^{K-5i^]r + 2(r.r20r2A 



1 f r r • r (r • r) r 



where the first term is the Keplerian acceleration, the second term has its origin 
in the perturbing potential of tidally deformed bodies. The term in the second 
row emerges due to the rotational deformation of both objects. The perturbing 
acceleration in the third row is for the relativistic contribution, while the last 
term comes from the dissipation of a mechanical energy. In spite of very different 



approaches used in (Eggleton et al 1998) and in this work, the final equations 



of motion, especially their dissipative part, are the same (see their Eq. 34 with 
f^F given by Eq. 45). According with the above notation Aq = mi .Rp ^l 0' ^i = 
TTlQ -Rf ^L 1 ^^'^ 
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4.2 The equations of the evolution of angular velocities 

The equations of the evolution of angular velocities of the star and the planet, 
denoted here by Qq and i?i are obtained through solving Eq. ( 11 ) with respect to 
fli {I = 0, 1). To the first order, one derives the following form of these equations: 



+ 



^1 
2 Ii 



ir-Ol)r+{r-Oi)r+{r-r)ni-5{r- Oi) (r ■ r) 



1 



{Ql X r) - r 



(50) 



where the first-row term is the same as would be derived directly from the Euler 
equations. The second-row terms emerge due to dependency of the potential energy 
on Ql . We note that the non-rigid-body contribution to the equations are the same 



order as the rigid-body term. The dissipative term stands in the last row of (501 
and it is exactly the same as in (Eggl eton et al|1998 l. To compare the results, see 
their Eq. 36 with fpp given by Eq. 45. 



4.3 Effects not included in the model 



The equations derived so far do not take into account a few physical and perturbing 
phenomena, which we would like to list here explicitly. First of all, we assumed that 
the star does not evolve in time keeping its internal structure constant. However, 
if to consider the evolutionary time scale, the stellar radius as well as the internal 
structure change: the time evolution of Rq and the mass distribution alter the 



inertia moment, which then may provide an additional term in Eq. (501. Hence, also 



the Love number is not constant over time. Evolution of the internal structure, for 
instance, the size of the convective zone, implies that none of coefHcients Aq remains 
constant. Also a significant stellar wind would be responsible for yet additional 
term in Eq. ([SO]), (e.g., Dobbs-Dixon et al„2004{ [Barker and OgiIvie|p009J . 

The same discussion applies to the planet. Its radius and internal structure also 
may vary over time, which is mainly attributed to the tidal and thermal influence 
by the parent star (e.g.,|Gu and Ogilvie|2009| [Miller et al|2009| [Hansen|20lol [Arras 



and Socrates[[2010 Ibgui et al|[2011 ). However, the tidal evolution of the interiors 
of hot-Jupiters is not well recognized. Moreover, very close-in planets revolving 
around their stars with periods of the order of one day, may be close enough to 
loose their mass through the Lagrangian point L3 ( Li et al|2010 l, as they may fill 
up the Roche lobe. 

In the considered model, we assume, that the material is purely viscous, which 
means that there is no unelasticity nor rigidity in the bodies. The model ignores 
the so-called yl-effect, which leads to differential rotation in the stars (see, e.g.. 



Kitchatinov and Rudiger 1993 Kueker et al 19931. However, this effect is well 



studied for single stars, it is not understood yet how a close planetary companion 
can tidally affect the differential rotation of the star. And on the other hand, the 
tidal influence of the star on the planet's differential rotation is not studied well 
enough. Here, we omit this effect. 
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As it was noted already, the dissipation parameters cr; (and also A/) depend on 
the tidal frequency, ui^. After (Eggleton et al 19981, we assume, that the energy 



is dissipated due to turbulent convection in the objects (thus, we consider the so- 
called equilibrium tide model). Therefore, u; parameter depends on the value of the 
turbulent viscosity coefficient ut (see Eggleton et al||l998 Eq. 113). We assumed, 
that ft is constant, i.e., does not depend on wt, which is not true. For small uj^, 
i/^ is nearly constant, while for larger tidal frequency, it is usually believed to 



be proportional to the inverse of its square, (see, e.g., Goodman and Oh 1997 



Terquem et al 1998). Moreover, ut is not isotropic, and its value depends on the 



Coriolis number (Kueker et al 19941. While for the star of the known internal 



structure, it is possible to calculate u^, for the planet it is usually not possible due 
to very poor knowledge of it's interior. In our simplified model, we do not take 
this effect into account. 

The equilibrium tide is not the only contribution to the energy loss. The other 
mechanism relies on the damping of the internal waves excited by tidal force 
(the so-called dynamical tide). The waves may be damped by radiative diffusion, 
convective viscosity as well as non-linear breaking (see, e.g.. Barker and Ogilvie 
2010 and references therein). This mechanism produces different dependence of 



the dissipative parameter, then the equilibrium tide. In the literature, the so-called 
quality factor Q (or rather the modified quantity Q') is used to express the energy 
dissipation rate. If we find formal relation between a and Q, the equilibrium tide 
model leads to Q (x uui. Dynamical tide model, on the other hand, gives the inverse 
relation. The dependence is in fact much more complex and is a subject of many 
studies in the literature (Zahn||1975 Goodman and Dickson||1998 Goodman and 



Lackner||2009[ |Ogilvie||2009[ [Barker and Ogilvie||2010[ |Efroimsky||20n| 



5 Averaging over the mean anomaly 



The orbital equations of motion have the following general form: 



-^r + f. 



(51) 



where the first term is for the Keplerian acceleration and f is a perturbation of 
small magnitude as compared to the Keplerian term. Therefore, the planet moves 
on weakly perturbed elliptic orbit. The time scales of variation of its elements 
are a few orders of magnitude longer than the orbital period. Hence, one might 
perform the averaging of the equations of motion over fast evolution (i.e., over the 



mean anomaly) making use of the averaging principle [see, e.g., Arnold et al ( 1993 1, 
( 1995 1, §51, §52]. Nevertheless, at first we should verify whether the 



Arnold 



§6.1 or 

evolution of Qq and Qi occur in slow-enough time scale. Because the magnitude of 
the moment of inertia of the planet is much smaller than that one of the star, i.e., 
/i <C /q, the planetary angular velocity vector v arie s faster than the stellar one. 
With the help of a relatively simple analysis of (|50[p| one finds that the relative 
variation of f2i during one orbital period has its order of: 



\Oi\ 

~n7 



p. 



orb 



5 7r mo j Ri 
2 mi \ r 



n 



^ The analysis, i.e., an estimation of the order of magnitude of \Q\ \ , is done by setting e = 0, 
the inclination between Q\ and r to j and physical parameters to their typical values. 
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where n is the mean motion. For a one-day orbit and fast rotating Jupiter-hke 
planet (i7i ~ n) , the above quantity is of the order of 6 x 10^^. For wider orbits, 
it obviously decreases. Thus, we have shown that the assumptions of the averaging 
principle are fulfilled at acceptable level even for rather "extreme" systems. Never- 
theless, we should keep in mind that this level of approximation is only acceptable 
when one would like to compare the results of the mean and exact model. 



Having the equations of motion in general form of (511, one may find the 
following equations for the evolution of the orbital angular momentum and Laplace 
vectors, respectively h and e: 

h = rxf, e= -|f X h + r X (r X f)|. (52) 

Keplerian orbital elements, such as the semi-major axis, eccentricity, inclination, 
longitude of ascending node and argument of pericenter, may be easily derived 
from the components of h and e. We choose this representation instead of the 
classical Gauss planetary equations, because it has no singularity with respect to 
the inclination. 

If the motion is considered to as the Keplerian one, all angles h, e, Qq and 
fii are constant and the fast vectors r and r may be expressed as the following 
combinations of e and h x e: 



_ e hxe ._. e . hxe , , 

r = XQj-b " + J/orb 7 5 ^ — ^orb ~ + 2/orb r > (5'jJ 

where ^orb ^'^'-^ 2^orb ^^^ Cartesian (a;, y)— coordinates in the Kepler frame, i.e, 
^orb ~ ^ COS!/ and j/orb ~ *" sin I/, where v is the true anomaly. Thus, the right- 
hand sides of the equations for h, e, Oq and Oi are some functions of ly, which are 
parametrised by constant vectors h, e, Qq and Hi. 

According with the principle of averaging^ the equations of motion are aver- 
aged out over the mean anomaly M, under assumption that the orbital elements 
are fixed over the averaging period. Formally, the procedure relies on calculating 
following integral: 



1 /■^'^ 1 /■^'^ (1-e^)- 



1 f 1 f 



(1 -|- e COS!/ 



,2 ' 



(54) 



where J- is averaged function. In the above formulae, we did a change of the 
integration variable, from the mean anomaly to the true anomaly. Note that we 



know ^ as a function of v rather then M (see the discussion in Migaszewski and 



Gozdziewski (2008)) 



^ We consider the problem using the averaging principle, although in general this principle 
is incorrect. Prom a technical point of view and when we don't need to know the relation 
between mean and actual elements, this principle is equivalent to the first order per turbation 
theory (for an overview on the classical perturbation theories see |Ferraz-Mello|2007|l . 



The generalized non-conservative model of a 1-planet system - revisited 



15 



By calculating such integrals over the right-hand sides of the equations of 
motion, one finds the secular equations that have the following form: 



^ = -2^ E^'-^i. 



2/3 



/=oa 



_9_ 

4/3 



^ aiAfj^2,i 



(55) 



;=o,i 



1=0,1 1=0,1 1=0,1 

Oi = ^AiJ',^i + ^aiAfj'2.i, 1 = 0,1- (57) 

Formally, variable names of h, e, etc., should be encompassed by square brackets 
{...), that denote the averaged values. However, it will be clear from the context, 
that after the averaging all functions should be understood as the secular or the 
mean quantities. The functions ^ are given by the following formulae: 



•^1,: 



J-. 



2,1 



-3,1 






|/i^ (e ■ Oi) {e X Hi) + [{e X h) ■ Qi] (e x h) x /?, |, 



3 



+h- 



' £5[{e X h) ■ Hi] {e X h)Y 



Ue^ {2fe^[(e.O,)(exh)-^ 



i e 



+ [2/i^e^r2f -7/i^(e-/2;)^-5[(ex h)-0,]^l (ex h) 
+2e^/i^[(ex h) ■«,]/?,+ 4 /i^ e^ (e ■/?;) (^/ x h) |, 



,3 /'i „2\ 



■ f 5 (e X h) , 



■7^5 = \_^ ' (e X h) , 



■6,i 



/lie 



|l8 ^^ fie + ft^ fs [ (e ■ fi,) h - 11 (h • /?,) e] |, 



where £i = fi(e) are the following functions of eccentricity 



o 1 15 2 15 4 5 6 



15 2 45 4 5 6 



^3 = ^+r' + r'' 



£4 = l + 3e +„£, 



c 1 I 3 2 1 4 

^^ = ^+2" +8"- 



All these functions are equal to 1 for circular orbits and are greater than 1 for 
elliptic orbits. To verify the correct form of the averaged equations, we compare 



the time-evolution of the unaveraged system, Eqs. ( 49 1 and ( 50 1 with the evolution 
of the secular system, Eq. (55l-(57l. Figure IT] illustrates the results of that com- 
parison as plots of the azimuthal angle of Q\ in the orbital reference frame, <jf>i(i) 
(the left-hand panel), and the argument of pericenter uj(t) (the right-hand panel). 
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Fig. 1 Temporal evolution of angles 4)-i(t) (the left panel) and u){t) (the right panel) derived 
in terms of the solution to the full equations of motion (black dots) and to the equations 
of the first- averaged system (gray curves). Main parameters of the system are mo = Ituq, 
mi = Imj, Ro = IRq, Ri = lRj, a = 0.02 au, e = 0.1, Tj-^t o = ^d . ^l-ot 1 = 1 d. 




15 t [yr] 20 



Fig. 2 Evolution of Q\/n(t) (the left panel) and 9\{t) (the right panel) derived by a solution 
to the full equations of motion (black dots) and the equations of the first averaged system 

"^ and 



(gray curves). Parameters of the system are mo = liriQ, ra\ = lm,j, Rq = IRq, Ri 



■ 0.02 au, e = 0.1, T'j.q^ q = 5 d, Tj.q)- j^ = 1 d. The dissipative coefficients are Aq 



10 



Ai = 10 . Red line is for quasi- equilibrium value of the Oi/n ratio. It is explained further in 
the text. See subsection 7.4 and Eq. (89) for details. 



The integration of both the non-averaged and the averaged equations of motion 



were done with the help of the Taylor integrator (software package by Jorba and 
Zou||2004l. As we can see, the results agree very well, although the precession of 



the planetary spin is relatively fast, with the period of about 100 days. We will 
discuss further in this paper, that this time-scale is the fastest one after the or- 
bital period, and the next one, in terms of the magnitude, is the characteristic 
time-scale of the pericenter advance. The experiment described above concerns 
the conservative system, with Aq = Ai = 0. The next figure [2] shows the results for 
the non-conservative models, with Aq = 10~^ and Ai = 10 (the parameters are 
unrealistically large, just to shorten the computations time). The left-hand panel 
is for the evolution of the Hi/n ratio. As we will show later on, for an eccentric 
orbit, the rotational frequency does not converge exactly to the mean motion n, 
but rather to some larger value. The left-hand panel illustrates the time evolution 
of di (angle between Oi and h). Clearly, it decreases to zero. Both these quanti- 
ties evolve in relatively short time-scale; let us recall that Ao,Ai parameters are 
by 3 — 4 orders of magnitude too large in this test as compared to likely values of 
the real systems. Again, the numerical test shows a perfect agreement between the 



The generalized non-conservative model of a 1-planet system 



revisited 



17 




0.08 



0.04 



9 t [kyr] 12 



Fig. 3 Evolution of e(t) (the left panel) and a(t) (the right panel) derived by solving full 
equations of motion (black dots) and the equations of first-averaged system (gray curves). 
Parameters of the system are mg = Ituq, mi = l"ij, Ro = 1^0) ^1 = l^Ji '^ = 0.02 au , 
e = 0.1, Tj.Q^- Q = 5d, Tj.q)- ]^ = 1 d. Dissipative coefficients are Aq = 10~^ and Ai = 10~^. 



results derived from both models. Eccentricity e and the semi-major axis a evolve 
in longer time-scale. Figure Is] shows plots of e(i) (the left-hand panel) and a{t) (the 
right-hand panel) for the same parameters as taken in the previous experiment. 
Again, the agreement between the exact and the mean models are excellent. 

As we have seen on Fig. [l] the planetary spin evolves in the time-scale that 
is only two orders of magnitude longer then the orbital period. It is still short, as 
compared to the full time-scale counted in Gyrs. The integration step would be 
then still too short to study the long term dynamics of the system CPU-effectively. 
In the next section, we will show that angle (j>i may be treated as the second fast 
angle, and it is possible to perform the second averaging of the secular system. 



6 Time-scale of the evolution and the second averaging 

Our next goal is to estimate the time-scale of the evolution of the conservative 
model. In this section we fix ao = cri = 0. That implies the following properties of 
the equations of motion: 



h-h = 0, e-e = 0, /2o ■ ^o = 0, /2i ■ /2i = 0. 



(58) 



It means that the conservative system possesses at least four first integrals e, h, Oq, f^i- 
The constant eccentricity and the magnitude of orbital angular momentum imply 
also the semi-major axis (a) constant. Moreover, it may be easily shown, that the 
total angular momentum L = /3h -|- /pi^o + ^i^i is a constant vector. Another 
constant quantity is e ■ h = 0. Therefore, the initial set of 12 scalar equations 
of motion may be reduced with the help of the 8 first integrals, so it should be 
possible to transform it to the set of four scalar equations. 

Nevertheless, we do not attempt to write them down, but rather to make use 
of the properties of the conservative secular system, in order to simplify the non- 
conservative system. At first, let us make a simple estimation of the time scale of 
the evolution of vectors h, e, Oq, ^i in the conservative system. We choose the 
following (say, representative) values of its parameters: mo = IniQ, mi = Imj, 
a = 0.02au, Ro = IRq, Ri = l/?j, ^o = 27r/(5d), ^i = 27r/(ld), A;l.o = 0.03, 
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kj^ I = 0.3, kq = 0.2, Ki = 0.5. We also take into account the dominant term only. 
We obtain: 






Por^---^(-] (^) 77^^^-8x10- 1-e^ \ (59) 



— -Porb ~ ^O TT 



e 



mi\a J (l-e2)5 V ^ ' ^ ^ 

72^^-b~^- — -^-j (1-ej ~3xl0 (1-ej ,(61) 

Jii 2 mi Ki n \ a J \ I V / 

As we have shown, Q\ becomes the fast vector after the first averaging. Then, we 
may express this vector in the Keplerian frame, assuming that h, e are constant 
vectors. We obtain: 

/ G h X G h \ 

Q\= Q\\ co%d>\ sin 6*1 — hsinoii sin^i hcos^i — I , (63) 

\ e he hj 

where 9i is an angle between h and Q\, and 01 is azimuthal angle in the orbital 
plane measured from the direction of vector g. Using Eq (571 for / = 1, one may 
find: 

01 = -^ ^ ^g ' ni cos6i, ^1=0. (64) 

The characteristic time-scale for 0i may be estimated as follows: 

r, -^"-^^OOdr " X^lEj^ mi Imprrot.l (1-e^)' 
■^1-0, U.02au>' V ^1 ; Imj mo Id cos^i ^^ 

The precession rate is constant in the conservative model. Nevertheless, for 9i f» 
7r/2, 01 may not change fast enough to be considered as the fast angle. That 
introduces a limitation of the analysis. To conclude, we stress that in some cases 
angle 0i may be not slow enough to make the averaging over the mean anomaly 
valid (that corresponds to too large, too close, too fast rotating planet), or it may 
be not fast enough to make it possible to perform the second averaging over itself 
(Uranus- type orientation of the spin vector). As we will show further on, the last 
limitation is weaker than the first one. 

After the second averaging, we obtain the following equations of motion: 

1 - ^ 9 / ,2-,__ , _. ,2, 



^^AoTuo-^[<yoAt,J'2,o + <^iAig2,ij, (66) 

e = --^ {A0F3.0 + A.Qs.i) - 15 f ^ + ^) Ta-^Fb 
Ap \mo miy c^ 

-7-5 f 0-0^0-^6,0 +0"1 ^106,1 J , (67) 

1 9 2 

^0 = 7ri^-4oJ^i,o + -pr croylo-^2,0, (68) 

2 Jo 4^0 
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Fig. 4 Evolution of f2(i/n{t) (the left-hand panel) and 8o{t) (the right-hand panel) derived 
from the equations motion of the first-averaged system (black dots) and the equations of 
the second-averaged system (gray curves). Parameters of the configuration are mo = 1"^0, 



mi = Imj, Rfi = IRq, Ri 
Dissipative parameters are Aq = 



= IRj, a = 0.02 au, e 
10-2 and Ai = lO^^, 



0.1, T, 



rot.O 



5d, T,ot,l = Id. 



where 






'3,1 



.6/1 oS'ii 



nf 5 cos^ 6)1-3 (e X h) 



06,1= ^^le ^ |18m gi-llfe^J^icosgigsje. 



(69) 
(70) 
(71) 



After the second averaging, the fastest angle is the argument of pericenter. Thus 
the time step-size of the integration increases significantly. To verify the correctness 
of this step, we perform a new experiments, which results are illustrated in Fig. |4j 
The parameters of the system are the same as in the previous tests. At this time, 
the plots show the evolution of i7o/n (the left-hand panel) and Gq (the right-hand 
panel). The agreement between the first averaged (black dots) and the second- 
averaged model (grey curve) are basically perfect. 



7 The third averaging and the dissipative equations of motion 



From equations ( 66 1-( 68 ) one can obtain equations governing evolution of a, e, Oq , 



and also i?i, 6i. They read as follows: 



9 
9 



or Y] <^l Af £q -£2(1- e^ ) 



4/?„8(i_,2)¥^f-^ 



Ts^'^l^ 



—^ COS 61 

n 



ISfi-llfs 



('-') 



5 n. 



cost/; 



(72) 
(73) 
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«0 



9 



o-Q Ann 



4/oa6(l-e2)^ 
9 aoAl 



2£:2cos6»o- (1-e^j' — (2£:4 - Psin6io) 



(74) 



47oa6(l-e2)'^ «o 



2 £2 sin 6*0 — 



f2o 



l-e' 



Pcoseo- i|2£:4sin6'o-p| 



Qi 



9 



(71 >lf n 



4/ia6(i_e2)6 
9 



2 f 2 cos 6*1 — f 4 ( 1 — e 



aiA? 



4/ia6(i_e2)b r2i 



smfi 



2£:2-£:4 1 



„-.})] 


' 


(75) 


^(l + cos^^i)' 


, (76) 


2)'^cosei' 
/ n 




(77) 



where 



Se 



V : 



^ ^ 31 2 _^ 255 4 _^ 185 6 _^ 25 g 
1 + T" +^^ +1^" +64^' 

^^ (e.r?o)" ^ ^^ [(exh).f?o]" 
e2 J72 sin 6*0 h'^ e^ i?2 gjj-^ g^ 



(78) 



These equations do not depend on the longitude of the ascending node nor on the 
azimuthal angle of Qq . 

A simple estimation of the time-scale of the equations of motion derived in the 
previous section shows that, after the second averaging, the argument of pericenter 
cj becomes the fastest angle in the system. It varies typically in a time-scale of 10^ 
years. To be more specific, and to find limitations to the choice of lo as the fast 
angle, we obtain a direct form of cj in the second averaged system (the conservative 
system). It is a sum of five terms, i.e., uj = cjs.o + cJ3_i -|- u)4,fi + oja.i -\- tJs, where 
particular terms of the sum correspond to functions ^3,0, -^3.1, -^4 (these are 
contributions form the star and the planet), ^5, respectively. We obtain: 



W3,0 



W3,l 



1^4,0 



W4,l 



W5 



Ao«o 



4/?a5^2 (l-e2)^ 

Ainl 

4/3a5^3 (l-e2)^ 



2 a cos ^0 + 5 cos ^0 ~ 1 



5 cos ^1 — 3 



a'5^(l-e2)- 

1 
15 /i 2 Ef, 



mo 
Ai 



a^(l^ 

3 
3p2 



e^y mi 



c2 02 (1 — e2) 
We introduce a new quantity: 



(79) 
(80) 
(81) 
(82) 

(83) 



I3h 

IqQq 



which is the ratio of the magnitude of the orbital angular momentum to the 
magnitude of the stellar rotational angular momentum. For typical parameters 
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Fig. 5 Critical inclination 9q as a function of a = q/(1 + a). Shaded areas correspond 

to retrograde rotation of the pericenter, while the white colored regions are for the prograde 

rotation of the pericenter. 



considered in this work, a is of the order of unity. Contributions form the tidal 
deformation of the objects iJJ4,o, i^4,i, as well as relativistic term ijj^ are always pos- 
itive. On contrary, the first two terms, stemming form the rotational deformation 
of the star and the planet, have some critical values of Jmut at which the sign of 
the rate frequency changes its sign (from positive to negative). 



7.1 The critical inclinations 

We may now introduce the critical inclinations 6q ' and 6[^^^ . The first quantity 
is a solution to the equation 0^3,0(^0) = 0, the second one is a solution to the 
equation dj3^i{9i) = 0. One can easily obtain the following expressions: 



^(crit) = - 1 (a ± Ta^Ts) , COS ef "*) = ± ^ ^ 



(84) 



The critical inclination 9^ ' has the following limits 



lim cos 61, 

a-»-0 



(crit) 



Vi 



lim cos 6, 



(crit) 



(85) 



Figure N shows graphs of cosSq " as a function of a = a/(l + a). In this 

paper we consider systems with a ~ 1 (or a ~ 0.5), thus magnitudes of the stellar 

^crit^ 
and orbital angular momenta are of the same order. Note that for &o ~ 6q , 

a contribution to a;3,o emerging due to the rotational deformation of the star is 

very small. If it was the only contribution, uj could not be considered as the fast 

ferity 
angle. Similarly, ii 61 ^ 6\ , the contribution of the rotational deformation of 

the planet is negligible. 
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Fig. 6 A time-scale of the rotation of periapses. 
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7.2 The time-scale of the rotation of pericenter 

To justify id as the fast angle in the second-averaged system, we estimate the 
time-scale t^i = 2 7r/tJ3 ;, r4 ; = 2n/ib4i { I = 0,1), ts = l-n jC:^. We obtain: 



T3,o « 10 yr 



mo 
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1^0 
mo 



1-e^ 



(86) 



where C; = cosS/ and rrotOi^lot 1 denote the rotation periods of the star and 
the planet, respectively. The typical (characteristic) time-scales were calculated 
for ^L = 0.03 and k-^ ^ = 0.3. 

The effect of the tidal deformation of the planet dominates over remaining 
perturbations in the range of the typical parameters. For a wider orbit, t^^i in- 
creases faster then rs.o or rs.i. Nevertheless, for a Sun-like star and a Jupiter-like 
planet, the effect of the rotational deformation of the bodies may be never the 
most important contribution to io. For a wider orbit, the dominating perturbation 
is the relativistic term. We have shown, that in the interesting range of physical 
and orbital parameters of the system, the argument of pericenter always increases 
(i.e., a; > 0), and may be treated as the fast angle in the third averaging. 

Figure |6] illustrates graphs of characteristic time-scale of the evolution of the 
conservative system. Thin, black curves are for the time-scale r^i, T41 and rs 
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discussed above. The dashed, red curve is for the join effect of two dominating 
perturbations, i.e., the general relativity and the tidal deformation of the planet. 
The time-scale due to the rotational deformation of the object are always longer 
than the former one, thus iu is always positive. This figure ensures us, that the third 
averaging (over oj), is valid in the considered range of parameters. Nevertheless, 
for fast rotating (periods of a one day) and large stars (with Ro ~ 2Rq), the 
time-scale ra.o rnay be shorter for some orbits (e.g., a = 0.1 an) than r4_i and T5. 



7.3 The time-scale of the evolution of the system 

Because of the small magnitude of the moment of inertia of the planet, it has to 
be verified whether the dissipative time-scale for Qi and 9i may be of the same 
order, as the time-scale of the rotation of periastron. For the completeness of the 
analysis, we estimate the time-scales of all variables, i.e., o, e, no,do, !^i,di- Because 
the right-hand sides of the equations of motion have rather complex form, we limit 
this estimation to the case of small e, and Oi/n ^ 1. For the typical parameters 
of the system, one finds the following characteristic time-scales: 



Ta,p « 5.5 X 10 yr 
Te,* ~ 1.5 X lO^yr 



10^6 / a \S flmeY ( mi y / 1 7?j 



Ai V 0.02 an/ y mo J \lmj J \ Ri 
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Ao V0.02au/ [imQj mi \ Ro 



10 
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lO-'^ f a yflmsY f mi Y f 1 Rj 



ra, « re„ ~ 1.2 X lO^r 
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mo / Y 1 mj / Y -Ri 
10'^ ( a y-i mo /lrnj\^ / li?© \ ^ 5d 



Ao V 0.02 an 7 1 m© \ m,i J \ Rq 
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where for any quantity denoted as X, the time-scale is defined as tx = X/\X\. 
The equations of motion for a and e are sums of terms representing contributions 
due to the energy dissipation in the star and in the planet, i.e, a = a^^' + a^'^\ 
e = eS*' -t-e'-P'. We calculated these time-scales for contributions coming from each 
body separately. The dissipative evolution of the planetary spin occurs relatively 
fast. However, one should keep in mind that these estimates were done under the 
assumption of Qi/n <C 1. This means that we neglect terms with Qijn in the 
equations of motion. This is not correct in general, and that simplification has 
been used only to derive the time-scales. Moreover, the magnitude of Ai is not 
well known. A concise discussion of the time-scales present in the system may 
be summarized graphically in Figure [Tl It illustrates a separation of particular 
time-scales for Keplerian motion, precession of the planetary spin (or frequency 
associated with the angle (/)i), rotation of the pericenter (a frequency of variations 
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Fig. 7 Time-scales of the Keplerian motion, precession of the planetary spin, rotation of the 
pericenter, and the dissipative evolution of the planetary spin. Parameters have typical values, 
'^rot,! = 1 d. 



of Lj), and finally, the time-scale of the dissipative evolution of planetary spin. 
According with formulae derived for r^^j , we included also the dependency of the 
mean motion, see Eq. (761 with e = 0,9 = 0. For extreme systems, with a one- 



day orbital period, the mean motion may be comparable with the period of the 
rotation of the pericenter, and Oi,9i may be considered as the fast quantities, 
similarly to uj. 

Nevertheless, the equations for Hi and 9i do not depend on lj, and the equa- 
tions for Oq and 9o, which depend on uu, do not depend of i7i, 6i. Therefore, we can 
average out the equations for the evolution of the stellar spin, without considering 
the evolution of the planetary spin, even if it varies in comparable time-scale. 



7.4 The third averaging and the quasi-synchronization of the planet's spin 



As we have shown, the time derivatives of Hq and 9q depend on the fast vector 
e, i.e., they depend on the fast angle uj. One can easily find that (V) = £4sin^o, 
thus, after the third averaging, the equations for Hq and ^o read as follows: 



^0 



aoAQn 



'iIoa<^{l~e^y 



2£2 cos 6)0 - £4 ( 1 - e^V — fl -h cos^ 6»o^ 



(87) 



00 Al 
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sin 6*0 



2£2-£A[l-e^y —\ coseo 



The spin of the planet evolves much faster than the other quantities, i.e., a, e, Qa,9Q- 
Assuming that these parameters are constant, one may find that Q\ and 9i tend 
to an equilibrium. There exists only one such an equilibrium, which is linearly 
stable, i.e.. 



fii 



> n, 
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£4(l-e2)^ 



h\ =0. 

eq 



(89) 
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Thus, the rotational velocity of the planet tends to n only for circular orbits, while 
for e > 0, the equilibrium value is larger than the mean motion. That is why we 
call this state as quasi- synchronous, rather than the synchronous one. For e > 0, 
the energy is still dissipated in the planetary interior, even if the spin has reached 
the equilibrium orientation. 

Time-scales Ta^p, Te.p are relatively short, ~ 10^ years for a one-day orbit. Yet 
they were obtained under the assumption that the planetary rotational velocity 
remains separated from the quasi-synchronous state, which is true only for t < tq^ . 
When the planet is located at this state, its contribution to a and e may be smaller. 
Putting the equilibrium values for fi\ and Q\ to Eq. ( 72 1 and ( 73 1 , we find the 
following expressions: 
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(90) 



(91) 



where 



P 1 ^ 45 2 ^ s 4 ^ 685 6 ^ 255 8 ^ 25 



10 



However, the tides in the planet cause a fast orbital decay; it acts in relatively 
short time-scale of the order of re,p. After circularization of the orbit, the tides in 
the planet do not contribute to a. We notice here, that in multi-planet systems the 
eccentricity as well as inclination of the orbit vary in the conservative time-scale 
and the spin of the planet may unlikely reach the equilibrium. Therefore, only in 
single-planet systems, it is possible to use the simplified equations of motion (i.e.. 



the equations with fixed 



and Q\ 



a 



lleqj 



8 The dissipative evolution of the system 



The final set of the equations of motion, i.e., Eq. (87 1, (881, (90), (911 has an 



integral of the magnitude of the total angular momentum, L = |L|, where 



L = /3h + /o/2o 



(92) 



To describe the evolution of the system, one has to solve three equations of motion, 
i.e., a, e and ^Oj for a- fixed value of L. Still, we cannot write down the solution 
to these equations, and we have to integrate them numerically. Moreover, we may 
study some particular solutions, like the equilibria, analytically. 
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8.1 The stability of the equihbrium 

The final approximation of the evolutionary equations of the system possesses one 
type of equilibrium. It corresponds to a state, in which the mechanical energy is 
not lost, i.e., e = 0, Hq = n* = Y/i/ai?,6'o = 0, where a* (or n*) may be treated 
as a parameter of a family of such equilibria. For some a*, the equilibrium may 
be stable, while for some other value - unstable. To study the stability, we write 
down the linear variational equations in the vicinity of the stationary solution as 
follows: 

. _ 27aoAl 9aoAl 
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are 


the variations of a, e, i7o, 6 



(93) 

?o, respectively. Equations 
for 5e and Sg^ are separable, and they admit simple solutions, i.e., both these 
quantities decrease exponentially to 0: 

5e(i) = 5e(0)exp(Aei), 5g^(t) = 5g^{Q) e^p{\e„t), Ae<0, Ae„<0. (94) 

The equations for 5a and 5q^ are conjugate each to other. To solve them, we define 
a new quantity x = Oi^/n and write down the relevant variational equation: 

QctqAI (Pal 



This equation also admits a simple solution: 5x{t) = 5x{Q) eyip{\xt), where Ax may 
be either positive or negative. The stability condition reads as follows: 

A. <0 ^ ^>3. (96) 

Jo 

For typical values of the physical parameters, the equilibrium is stable if 

The orbital period corresponds to a = 0.074 an is ~ 7.4 d. Thus, for a Jupiter- like 
planet orbiting a Sun-like star, with orbital period shorter then ~ 7.4 d, the planet 
does not tend to the equilibrium. When the equilibrium is unstable, the initial 
condition Qq < n forces the planet to fall down onto its host star. 

Considering the time-scale of the tidal evolution of orbits with semi-major 
axes larger than this critical value, one may conclude that the equilibrium may 
be unlikely attained by a Jovian planet. Moreover, in a more realistic model, this 



equilibrium does not exist at all (Barker and Ogilvie 2009), because additional 



effects, like the angular momentum dissipation due to stellar winds, may decrease 
the rotational velocity of the star. 
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Fig. 8 Evolution of Qa/n{t) (the left-hand panel) and a{t) (the right-hand panel) derived 
through the solution of the equations of motion of the second averaged system (black dots) and 
the equations of third-averaged system (gray curves) . Parameters of the system are mo = 1 "^0 , 
mi = liTij, Rq = IRq, Ri = IRj; o. = 0.02 au. Dissipative coefficients are Aq = 10~^ and 



Ai 



10" 



The initial system is close to the equilibrium state. 



Figure IS] illustrates the evolution of the system initially located in the vicinity 
of the equilibrium. Black dots are for the solutions to the equations of motion of 
the second averaged system, the grey curve is for the third-averaged system with 
the spin of the planet in the synchronous state. This test shows that the model is 
correct. It also illustrates the behavior of the system near the unstable equilibrium. 
The system with initially Hq > n excites both the semi- major axis and the Qq/u 
ratio. Still, the rates of variations of a and J7o/n do not suppress dissipation of the 
total energy. On the other hand, for initial i7o ^ n, the planetary destiny is to fall 
down onto the parent star. 



8.2 Parametric study of the dissipative evolution 



Here, we present the results of the analysis of the evolution of the system for 
various initial conditions. We adopt typical physical parameters of the system as 
follows. The masses are mo = \mQ and mi = Imj, for the star and the planet, 
respectively. Their radii are _Ro = 1 Rq and Ri = 1 -Rj . The mass distribution 
is described by polytropic models with indices no = 3 and n\ = 1.5. Parameters 
of the energy dissipation rates are Ao = 5 x 10^^ and Ai = 5 x 10~^ in the 
first case (results are presented on Figs. ^ and 10 1 and Ao = Ai = 5 x 10^^ in 



the second case (Figs. 11 and 12 1. We also consider a few sets of initial Tj-ot.o 
and 6q. Figure [9] illustrates the results derived for the initial rotational period of 
the star Tpot,o = 5 d, and for three initial values of the inclination between the 
stellar equator and the orbit (from the top to the bottom row, it is 0°,60°, 120°, 
respectively). Colors encode the final a (the left-hand column) and final Trot.o (the 
right-hand column). 

At first, lets us study panel (a), obtained for initial 6*0 = 0. For the initial 
(a, e) in the white region, the planet falls down onto the star during time which 
is shorter than 5 Gyr. The black, thick curve is for the border of survival. For the 
initial conditions located at the (a, e)-plane, on the right side of this curve, the 
planet does not fall onto the star during the whole time of integration 5 Gyr. Other 
black (thin) curves are for initial conditions for which the planet survives during 
10^ yr, 10® yr and 10^ yr, respectively. As we can see, for parameters adopted 
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Fig. 9 Maps of the final state of the system as a function of the initial semi-major axis and 
eccentricity. Panels in the right-hand column are for the final a, the right-hand column is for 
the final rotational period of the star. Parameters of the system are the following: mi = IrriQ, 
mi = 1 mj, Rq = 1 Rq, Ri = 1 Rj. Polytropic indices of the star and the planet are equal to 
3 and 1.5, respectively. The energy dissipation constants are Aq = 5 X 10~^,Ai = 5 X 10~^. 
The initial Tj-ot^ = 5d. The initial angle 6*0 is 0°,60°, 120° for the top, the middle and the 
bottom rows, respectively. The integration time is 5 Gyr. 



in this experiment, the border is placed at ~ 0.046 au for initially circular orbits 
and at ~ 0.068 au for initial e = 0.5. Clearly, planets with initially larger e migrate 
towards the star faster. As we already noticed, the contribution to a stemming from 
the energy dissipation in the interior of the planet is typically larger or similar to 
that one due to the dissipation in the star. Nevertheless, this term is proportional 
to e and vanishes for circular orbit. The eccentricity is damped in similar time- 
scale than a for corresponding values of Aq, Ai, thus the dissipation in the planet is 
important not only at the early stages of the dissipative evolution. The evolution 
of the eccentricity may be also read from the discussed figure. Contours plotted 
with yellow thin curves are for the levels of constant final e, which are 0.001, 0.1, 
0.2, 0.3 and 0.4, respectively. As we can see, for initial conditions located close 
to ~ 0.01 au on the right-hand side of "the survival curve", the eccentricity is 
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Fig. 10 The same as shown on Fig. 9 but for initial Tj-q^ q = 2 d. 



damped significantly, while for larger initial a, its variability remains small. For 
a > 0.074 au, the dissipation time-scales re are longer than the integration time, 
and the eccentricity is not altered significantly. Also the decay of the orbit is much 
slower. 

The second panel (b) shows the color-map of the final Tj-ot o- The black curves 
are again for the time of planetary survival. The border of 5 Gyr is also a border 
which divides the map of the final state of the star onto two parts. In cases when 
the planet falls down onto the star (on the left-hand side of the border curve), 
the rotation of the star becomes faster. For the initial conditions near the border, 
the rotational period decreases from 5 days to ~ 2.2 days. It may be understood, 
because the initial orbital angular momentum of the planet increases for larger 
initial a. For initial a > 0.074 au, for which the planet survives, and a variation of 
a is not large, also a change of Hq is not large. 

The next two rows show the results for initial Oq = 60° (the middle row) and 
120° (the bottom row). In the case of panels (c) and (d), the differences between 
these results and the previous ones are not significant. The border of planetary 
survival shifts slightly towards larger initial a. The minimal values of the final 



30 



Cezary Migaszewski 



0.02 




0.02 



0.02 



0.04 



0.06 



0.08 o|au| 0.1 



0.02 



0.04 



0.06 



0.08 o|au| 0.1 



Fig. 11 The same as in Fig. |9J but for Aq = 5 X 10 ^, Ai = 5 X 10 ^. The initial rotational 
period of the star is Tj.qj g = 5 d. 



2rot.o for ^0 = 60° increase also slightly, ~ 3 d near the border. This is an effect of 
addition of two vectors, i.e. /? h and Iq Oq, which are mutually inclined. The angular 
momentum exchange that may manifest itself in the decrease of Tj-ot^o is smaller 
here than in a coplanar system. The inclination 9o does not change significantly 
during the evolution. The yellow contours shown at panel (d) indicate constant 
levels of the final Oq = 50°, 57°, 59°. The systems near the survival border change 
their 9o by about of 10 degrees, while for systems located close to ~ 0.02 au, on the 
right-hand side of this border, the inclination is almost constant. The observed slow 
rate of inclination dumping agrees with the determined inclination in hot-Jupiter 
systems (e.g., [Triaud et al||2010| |Pont et all|2010| ). 

For the retrograde orbits, i.e., with initial Oq > ^, here 6o = 120°, the final state 
of the system is different. The survival border shifts again. But most significant 
differences may be visible in the final rotational period of the star. The picture is 
somehow a negative with respect to those obtained for Oq = 0° and 6o = 60° . The 
larger final Ti-ot.o ^^re for those initial conditions for which the planet falls down 
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Fig. 12 The same as in Fig. |9J but for Aq = 5 X 10 ^, Ai = 5 X 10 ^. The initial rotational 
period of the star is Tj-q^ g = 2 d. 



onto the star or has survived but stays close to the border. The difference may be 
explained due to the z -components of initial /3 h and Jo Oq have opposite signs. 

Again, the inclination of configurations which end with a > does not change 
significantly. The evolution of the eccentricity in all three cases looks like very 
similar. The main difference between the studied cases concerns the final Trot.o- 
The next figureflOlillustrates the results derived for parameters used in the previous 
experiment; only the initial rotational period of the star is now Ti-ot.o = 2d. The 
view of the final state of the system is similar to that one presented in Fig. |9J 
There is a difference, however, that relies in a shift of the survival border towards 
smaller initial semi-major axes for Oq < ^, and towards larger a for 6o > ^. It may 
be understood rather easily. For the faster rotation of the star, the border Qq = n 
shifts at the (a, e)-plane towards larger initial a. A contribution proportional to 
the ratio Qq/ti implies an increase of a, e, 6*0 and a decrease of Qq; it acts opposite 
to the remaining terms. Larger Oq/ti, for some selected initial condition, means 
that this term is more significant from these terms. It manifests itself also in the 
fact that for initially coplanar configurations, there exist regions in (a, e)-plane, in 
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Fig. 13 Time evolution of i7()/n(t) (tlie left-liand panel) and a{t) (the rigtit-tiand panel) of 
the third-averaged system. Parameters of the system are mo = 1 m©, mi = 1 mj, Ro = 1 Rq, 
Rl = IRj, a = O.OSau. Dissipative parameters are Aq = 10~^ and Ai = 10^*. The initial 
system is close to the equilibrium (i.e., e ft: 0, 6o Ri 0, i?o ~ !)• Grey curves are for the initial 
QqIu < 1, the black curves are for flo/n > 1. 



which the eccentricity is excited to larger values. For systems with 9o > -k 12 the 
survival border shifts towards smaller a. It is quite clear, because the term Hq/ti is 
multiplied by cos^o which is negative for the retrograde orbit. This term causes a 
decrease of both a and e, and has larger magnitude for greater Qq (shorter Tj-ot o)- 

In the next two figures |11| and |12| we present a similar study to the previous 
experiments, but derived for larger values of the energy dissipation constant Ai = 
5 X 10~^. The results are illustrated in the same manner. The only difference relies 
now in significantly smaller final eccentricity than it was derived in the previous 
experiments. Again, it may be explained if we recall that the contribution ep 
increases by two orders of magnitude with respect to the previous tests. For large 
Ai it dominates over e* and forces a fast dumping the eccentricity. On the other 
hand, the term ctp is important only at the early stages of the evolution, when 
e is significantly different from 0. When e becomes very small, the term ap does 
not accelerate the orbit decay. The borders of survival as well as maps of the final 
Tj-ot ^^^ very similar to those ones obtained for smaller Ai. 

It should be noted here, that in none cases illustrated on Fig. [9]|12[ the system 
ended in the synchronous state. For a < 0.074 an this equilibrium is unstable, while 
for a > 0.074 au the time-scales of the evolution are too long to fix the system in this 
state. Figure [13] illustrates the evolution of the system which is initially close to the 
equilibrium for a = O.OSau (the rotational period of the star Trot,o ^ [5.4, 9.4] d). 
Although the equilibrium is stable, the system tends towards this state very slowly; 
the time on the x-axis is counted in terayears\ Moreover, for initial i7o/n < 1 (the 
grey curves), the planet likely will fall down onto the star, even if at the beginning 
the system seems tend to the equilibrium. It may be explained relatively easily. If 
the initial fip/n- < 1, then for a > 0.074 au this ratio increases, tending to unity, 
and simultaneously a decrease. When a decreases below the limit of stability, the 
equilibrium becomes unstable and the system evolves outwards the equilibrium. 



9 Summary and conclusions 



In this work we revisited the problem of the generalized model of a single-planet 
system including the mechanical energy dissipation. We derived the equations of 
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motion from the very basic formulation of the Lagrange equations of the second 
kind. In that approach, the potential is a function of the angular velocities of the 
bodies. Our derivation is different from the standard approach used to derive the 
rotational equations of motion of the rigid body, in which the potential depends 
only on the Euler angles describing the orientation of the rigid body, and does not 
depend on their time derivatives. We obtained the equations possessing a different 
general form, but in the particular case considered here they remain very similar to 
the Eulerian equations. Moreover, as we show, an additional term appearing in the 
right-hand side of the equation for Qi does not introduce any secular contribution. 

The derived equations of motion were averaged out over time-scales corre- 
sponding to the conservative evolution of the system. We analyzed all characteris- 
tic time-scales, and we showed that they may be ordered in terms of a hierarchical 
set of variables, which then may be treated as fast variables in a recursive av- 
eraging process. The fastest component of the evolution is the Keplerian mean 
motion of the planet. After the averaging over the mean anomaly, we obtained the 
first-averaged system. This set of equations describes the evolution in which the 
fastest variable is the precessing angular momentum vector of the planet. Then 
the second averaging was performed over the azimuthal angle of i?i in the orbital 
reference frame. In this way we obtained the second averaged system. Finally, the 
third averaging was performed over the argument of pericenter. Thus, to eliminate 
secular variability that occurs in the conservative time-scale, and to obtain the 
dissipative equations of motion, we had to average out the system over three fast 
angles, i.e, A4,(f)i,uj. The precision of the averaging has been tested at all stages 
of the procedure, and we demonstrate that this approach leads to correct results. 

Next, we have shown that the dissipative evolution of the angular momentum 
of the planet occurs much faster than the orbital evolution as well as a variability 
of the stellar angular momentum. The inclination of i?i with respect to h decrease 
to in the time-scale only slightly longer than the characteristic time-scale of the 
pericenter rotation, for a close-in planet. Simultaneously, Qi tends to an equilib- 
rium value which is > n, where the frequencies are equal only for the circular orbit. 
The equilibrium is stable and we can fix di and Hi at their equilibrium values. In 
this way we derived the final set of equations governing the long term evolution of 



the system that admits energy dissipation, Eq. (87), (881, (90) and (91). 

The obtained equations of motion describe a dynamical system with one equi- 
librium corresponding to the circular, coplanar and synchronized orbit, which is 
unstable for orbits inside certain critical distance from the star, and which is sta- 
ble for orbits outside this border. For a Sun-Jupiter system that border is on 
a ~ 0.074 an. Studying the evolution of the system for generic values of physi- 
cal parameters, we have shown however, that this equilibrium is unlikely to be 
reached. As we demonstrated, the final state of the system is a complex function 
of the initial conditions and physical parameters of the model. The derived equa- 
tions make it possible to study this problem very effectively, both analytically and 
in terms of the CPU overhead, because all the evolution of dynamical variables 
occurring in the intermediate time-scale related to the conservative corrections 
have been averaged out. 
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Appendices 



A A proof that IS given by formulae (45| fulfills the condition ([2|) 



In the considered case, the left-hand side of condition ([2]l reads as follows: 



i = l 



dE 



dE 



dE 

1h 



dE 



P 



dqi dr ds dp 

It is quite easy to show, that the following equality occurs, when E is given by Eq. |45| 



dE ^ dE ^ dE ■ 

r ■ h Oo ■ \- ^1 ■ =2E. 

dr dOo dQi 

Now, to prove a consistency of the Lagrange equations with a dissipative term in the considered 
problem, we have to show that 



dE ^ dE 

S ■ = ill) ■ 

i9s dOd 

(and similarly for p and Qi terms). It is again quite elementary and may be checked by using 
the above formulae for dE/ds and the relation between Q and s (equation l4jl as well as the 
fact, that the following equality occurs for unit quaternion s (i.e., Sg + s^ = 1): 



dsQ 
ds 



s 

so' 



Q.E.D. 



B Obtaining the equation (m) 

The equation in the middle row of Eq. iIm may be transformed into equation ^ by using the 
fact that Qo = Qo{s, s). For the component x we have (index in Qq is omitted): 



dsx 



dn df 
d Sx d Q 



dj_ 
dsx 



dn df 
d Sx d Q 



where / is for C or E. For the time derivative appearing in the Lagrange equation we obtain: 



d / dC 
dt \dsx 



d / dn 

dt \d Sx 



dC 

dn 



dn d /dC 
d Sx dt \dO 



Similar expressions may be obtained for y and z components and then we write the following: 



/ d dOx d dQy d dn^\ / dC \ 



'dC 
dt \~dl 



dt d Sx dt d Sx dt d Sx 
d d ffx d d fly d d flz 



d d Qx d d fly d d fl^ 
\ dt d s, dt d s, dt d s, I 



dj^ 

d fly 
di 

\dfl, ) 



/dilx dSly c)il2\ / d t)L \ 



dftx dffy d fiz 



\T. 



ds, ds. 



dt d fix 

d d£ 

dt d fly 

d dC 
\ litdfl, J 
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f dQx 



d Hy d Q, 



dC 




1 dE 


dt d Sx 

d ddx 


dt d Sx 

d dOy 


dt d Sx 
d dff. 


2 9s 


dt dSy 

d dQx 


dt dsy 

d any 


dt dsy 

d dn^ 




\dt ds^ 


dt dsz 


dt dsz 



Collecting these terms according to middle row of Eq. 



/ 1 dE \ 

2 dQx 

1 dE 

2 dfly 
1 dE 

V 2 ~dQl ^ 
we obtain Eq, 



0. 
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